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Abstract 

The purpose of this paper is to present asymptotically stable open boundary con- 
ditions for the numerical approximation of the compressible Navier-Stokes equations 
in three spatial dimensions. The treatment uses the conservation form of the Navier- 
Stokes equations and utilizes linearization and localization at the boundaries based on 
these variables. 

The proposed boundary conditions are applied through a penalty procedure, thus 
ensuring correct behavior of the scheme as the Reynolds number tends to infinity. The 
versality of this method is demonstrated for the problem of a compressible flow past 
a circular cylinder. 
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1. Introduction. In the present paper, we discuss boundary conditions for dis- 
sipative, wave dominated problems, exemplified by Burgers equation and the three- 
dimensional, compressible Navier- Stokes equations given in conservation form. The 
emphasis is on deriving open boundary conditions ensuring the continuous problem 
to be well-posed and on devising semi-discrete schemes for imposing these conditions, 
which can be proven asymptotically stable. The boundary conditions and the semi- 
discrete scheme are valid even in the limit of infinite Reynolds number. 

When addressing exterior, wave- dominated, dissipative problems, one is often 
forced to introduce an artificial boundary for computational reasons. This introduces 
the well known problem of specifying appropriate boundary conditions at the artificial 
open boundary. For purely hyperbolic problems, it is well known that enforcing these 
boundary conditions through the characteristic variables leads to a stable approxi- 
mation. However, for dissipative wave problems the procedure is considerably more 
complicated. 

Naturally, we must require the boundary conditions to lead to a well-posed con- 
tinuous problem. For wave problems of dissipative type, the problem must, in order to 
be compatible with weak boundary layers, remain well-posed even in the limit where 
the dissipation vanishes and the problem becomes purely hyperbolic. In addition to 
this, we wish the discrete approximation of the problem to be asymptotically stable, 
and that the boundary conditions are easily implemented. 

For general non-linear problems the issues of well-posedness and asymptotic sta- 
bility are very complicated, and for most problems relatively little is known. However, 
as discussed by Kreiss and Lorenz [1], we may, for a large class of operators, simplify 
the problem significantly if the solutions are smooth. It was shown that in this case 
it is sufficient to consider the questions of well-posedness and asymptotic stability for 
the linearized, constant coefficient version of the full problem. 

The energy method is applied to the linearized, constant coefficient version of the 
continuous problem in order to obtain energy inequalities which bound the temporal 
growth of the solutions to the initial-boundary value problem. This technique allows 
for handling such complex problems as the Navier- Stokes equations and is in general 
applicable to symmetrizable problems [2], 

The usual way to enforce the boundary conditions in the numerical scheme, once 
their proper form for the continuous problem is known, is to solve the equation in the 
interior of the computational domain, and then enforce the boundary conditions at the 
boundary points. However, this approach does not take into account the fact that the 
equation should be obeyed arbitrarily close to the open boundary. To circumvent this 
problem, Funaro and Gottlieb [3, 4], and Carpenter et al. [5] developed the penalty 
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method which enforces the boundary conditions, as well as taking into account the 
equation at the boundary. They showed asymptotic stability for the scheme applied to 
scalar hyperbolic equations and systems of hyperbolic equations. Don and Gottlieb [6] 
recently showed how this idea can help in applying the Legendre collocation method 
on Chebyshev grids. 

The proofs presented in this paper are all done for semi-discrete schemes. The 
relation between the stability of the semi-discrete and the fully discrete scheme was 
recently discussed by Kreiss and Wu[7]. 

The issue of well posed boundary conditions for the compressible Navier-Stokes 
equations was previously considered by Gustafsson and Sundstrom [8], Oliger and 
Sundstrom [9], and Nordstrom [10]. They all used the energy method to derive bound- 
ary conditions for the linearized, constant coefficient Navier-Stokes equations in the 
primitive variable formulation. Dutt [11] introduced an entropy function, which al- 
lowed him to derive boundary conditions for the non-linear problem, ensuring that 
the solution remains bounded in an entropy norm. 

The remaining part of this paper is organized as follows. In Sec. 2 we review 
some well known results on Legendre polynomials and collocation methods. Section 
3 discusses Burgers equation, and boundary conditions ensuring well-posedness of 
the problem are derived. We continue by proposing an asymptotically stable penalty 
method through which the boundary conditions are enforced. This scheme ensures the 
correct behavior even in the limit where the problem becomes hyperbolic, and may in 
general be applied to any non-linear scalar equation. The penalty method for scalar 
hyperbolic, parabolic, and linear advection-diffusion equations is briefly discussed, 
and the proposed scheme is evaluated by numerical tests. The importance of properly 
choosing the penalty parameter is addressed in Sec. 4, where we discuss the effect of 
the penalty method on the CFL condition when using explicit Runge-Kutta methods 
for time-stepping linear problems. The results from the linear analysis are s! h! own 
to carry over to the non-line 

2. Legendre Polynomials and Collocation Methods. The schemes which 
we analyze in the present paper are all based on Legendre collocation methods. This 
choice is merely dictated by a wish to obtain analytical results, and the methods ex- 
tend trivially to other collocation methods and even to finite difference /finite element 
methods. 

The Legendre polynomial of order N is defined as 

1 

P "M = WNi.^ X ‘-V ' 

where \x\ < 1. We will in the following only consider collocation methods, where the 
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collocation points are given as the Legendre-Gauss-Lobatto points, being defined as 
the roots of the polynomial (1 - x 2 )P* N (x). There is no known explicit formula for 
these roots. 

Associated with the Gauss-Lobatto points is the quadrature formula, stating that 
if f(x) is a polynomial of degree 2 N — 1, then 

N 

= 

k = 0 

where x k are the Legendre-Gauss-Lobatto collocation points, and the Gauss-Lobatto 
weights, u are given as 


/'/«)« . 


( 2 ) 


k N + 1 PN{*k) PN-\{*k) 
2 

^0 = 


1 < k < N - 1 




N(N + 1) 

For further details on the properties of the Legendre polynomials, we refer to [12]. 

In a Legendre collocation method, the function, f(x) } is approximated by a grid 
function, fk — f{xk)> where the grid points are the Gauss-Lobatto collocation points. 
Thus, we construct a global Legendre interpolant, J/Vj to obtain the approximation to 
the function; 


N 


{lNf){x) = ^2fkh k {x) , 


1=0 


where the interpolating Legendre-Lagrange polynomials are given as 


h k (x) = 


(1 - x 2 )P^(x) 


N(N + 1) (i - Xk)PN(xk) 


We note that by construction 


(7yv/)(x/t) = fk ■ 

To seek equations for an approximate solution, (J/v/)(z), to a partial differential 
equation, we need to obtain values for the spatial derivatives at the collocation points. 
This is done by approximating the differential operator by a matrix operator, with 
the matrix entries given as 

J) ki - h\(xk) • 

For the explicit expression of the entries, we refer to [13, 14]. 
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3. Burgers Equation. In this section, we consider Burgers equation 


(3) 


dU TT dU d 2 U . . „ ^ n 

U ^ ^ 0 I® 1^1 <> o , 


dt dx dx 2 
where e > 0. The initial condition is given as 

U(x,0) = f(x) , 

with boundary conditions of the form 


(4) 

(5) 


aU(-l,t)-06^ 
7 U(1 > t)+6£ ^ 


= 0 


X— — \ 
= 0 


a?=l 


When addressing the issue it is, as discussed in the introduction, sufficient to consider 
the linearized, constant coefficient version of Burgers equation 


( 6 ) 


dU dU _ &HJ_ 
dt + dx 6 dx 2 


\x\ < 1 t > 0 . 


Here A = Uo is the uniform solution around which we have linearized. Equation (6) is 
also known as the linear advection-diffusion equation. 

The four real constants, a, /?, 7, and S , in the boundary conditions, Eq.(4)- 
(5), may not be chosen arbitrarily, since the resulting problem should be well-posed. 
Bounds yielding a sufficient condition for well-posedness are given in the following 
Lemma. 


Lemma 3.1. Equation (6), with boundary conditions given by Eq.(4)-(5) is well- 
posed if one of the following conditions holds 

(i) : 0 = 0 , 6 = 0. 

(ii) : 0^0 , 6 = 0 and (e — A) + 2a/ 0 > 0. 

(Hi): 0 = 0 , 6 ^ 0 and (e + A) + 2*y/6 > 0. 

(iv): 0^0 , 6 ^ 0 and 2{e - \)i/6 + 2(e + A)a//3 + 4(ay)/(06) > A 2 . 


Proof Construct the energy integral as 


\f t \\Uf = -HU,U x ) + eW,U xx )=\ 


-XU 2 + 2eUU x ] 1 ^ -e\\U x \\ 2 . 


Here we have introduced 

(U,V) = fuvdx , ( U,U) = \\U \\ 2 . 
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Following the similar analysis done in [15], we use the following estimate 


-e\\U x W 2 <- £ -[U(l)-U(-l)r . 


Applying this, the condition for well-posedness becomes 


\jt m2 -\l~ xu2 + 2eUUx 


-l 




Condition (i) implies that U( — 1) = {7(1) = 0 such that 

<0 . 


For condition (ii) we obtain U(l) = 0 and thus 
yielding the condition 


£ — A + 2— > 0 . 

Likewise, for condition (iii) we obtain 

55» lr l |,£ -K e + * + 2 «D 1,2(1)20 ' 

showing that this choice yields well-posedness. For condition (iv) we obtain the fol- 
lowing condition 

\j t \\n 2 < (e - A + U\- 1) + eU(-l)U(l) - \ (e + A + 2 |) U 2 ( 1) < 0 . 

This is obeyed if 

e 2 — ^£ — A + 2—^ ^£ + A + 2~^ ^ 0 j 

implying 

2(e - \)i/6 + 2(e + A)a//3 + 4(a7)/(/?£) > A 2 . □ 


3.1. The Semi-Discrete Scheme. Equation (3) will be solved using a Legendre 
collocation method where the collocation points are the Legendre-Gauss-Lobatto points. 
This involves finding an Nth degree polynomial, u(x } t ), satisfying 


du du d 2 u 
dt ^ U dx € dx 2 


at x = Xk , 


ke[i...N-i\ , 


( 7 ) 
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in the interior. The boundary points are given by boundary conditions of Robin type 


/ \ n & u 

au{x 0 ,t) - (3e — 
du 

ru(x N ,t) + 6e— 


= 9\{t) 

= 92{t) 


x N 


> 


) 


where g\(t) and ^ 2 (0 are the boundary conditions. The traditional method of imposing 
the boundary conditions is to solve Eq.(7) in the interior and enforce the boundary 
conditions at the boundary points only. However, this approach does not take into 
account the fact that the equation must be obeyed arbitrarily close to the boundary. 
In addition to this, it has proven difficult to implement Robin boundary conditions 
consistently for non-linear problems. To overcome these problems, we follow the line 
of thought initiated by Funaro and Gottlieb [3, 4] and propose a penalty method 
for Burgers equation at the Legendre-Gauss-Lobatto collocation points, x = Xk E 
[0, . . . , N], as 


( 8 ) 


du du 
dt U dx 


where 


d 2 u 
dx 2 


Tl <3 (®) 


/ > a du 

aix(z 0 ,f) ~ 


x 0 


- 9i(t) 


t 2 Q + {x) 


1 u{xn, 0 + Se 


du 

dx 


- 52(0 


X N 


( 9 ) 


Q (*) 






(1 + 
2 ^( 1 ) 


These two functions have the property of being zero at all Legendre-Gauss-Lobatto 
collocation points, except at the two endpoints of the domain. Although Q~ and Q+ 
here are defined as delta-functions at the boundary, we may equally well chose other 
definitions. As shown by Don and Gottlieb [6], this approach may also be applied for 
implementing Legendre methods on Chebyshev grids. 

We note here that the penalty method as given by Eq.(8) combines the boundary 
conditions and the governing equation into one equation. When using the penalty 
method, the boundary conditions are not exactly obeyed at the boundary. However, 
the method remains spectrally accurate, as we will soon illustrate. One may also 
observe that the scheme is equivalent to the traditional approach for T \ , r 2 approaching 
infinity. 

In order to obtain the energy inequality, we consider only homogeneous boundary 
conditions. As discussed previously in [1], this is no restriction, since we may always 
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introduce a variable transform such the boundary conditions become homogeneous. 
In the following Lemma we state the bounds on T\ and Ti ensuring that the linearized, 
constant coefficient version of Eq.(8) is asymptotically stable. 


Lemma 3.2. Assume u(x,t) exists and let r b and r+ b be defined as 

T~ b = ^-L[e + 2/c - 2\J k 2 + en - l/2eu;|A|] , 

r + 6 = + 2k + 2\J k 2 + £K - l/2ea;|A|] , 

where k = wa/b and 

2 

“ “ N(N + 1) ’ 

is the Legendre weight at the end-points. 

Then if 

t ~ p < rj < r + i/3 , 

r " 6 < r 2 < r + 6 , 

then the linearized, constant coefficient version of Eq. (8) is asymptotically stable and 
the solution is bounded as 

Id...., ^ fdu \ 2 

2 S wu < g <■* • 

Proof We start be defining the discrete, weighted scalar product as 

N 

{u,v) N = ^2u(x k )v(x k )u} k ,(u,u) N = ||u||n . 
k = 0 

and note that since we are using a Legendre collocation method, we have, through 
Eq.(l), the identity 

{u,v x ) N = (U,V X ) . 

This makes it straightforward to apply partial differentiation. Following the results 
stated previously, it is sufficient to obtain the energy estimate for the linearized, con- 
stant coefficient version of Eq.(8); 

“IMI n = ~^[ uu x\ X -\ + £ [uur]l-i - £ IK||w 

— T\uu{ — l)[cm( — 1) — /3eu x (-1)\ — r 2 t^n(l)[ 7 u(l) + Seu x ( 1)] . 
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Here the subscripts designate differentiation and cv is the Legendre weight at the 
endpoints (Eq.(2)). Using the quadrature rule allows for rewriting 

7V-1 

IMI N = + «x(l)w + Yj u U x k)vk . 

k = 1 

Contrary to the approach followed by Funaro and Gottlieb [3, 4], we recast the prob- 
lem of stability into an algebraic eigenvalue problem. For the present problem, this 
may seem an additional complication. However, we find that for more complicated 
problems, this approach greatly simplifies the proofs. 

Isolating the terms contributing to stability at each boundary, we obtain two 
conditions for asymptotic stability; 

< 0 , u+?f + u+ < 0 , 

where u_ = [u(-l), u x {-l)] T , u+ = [u(l), u x (l)] T and 

A — 2au/ri — e(l — Put\) 

- (3u)T\) - 2 eu) 

-A — 2'yUJT 2 £(1 - 6 UJT 2 ) 
e(l — 6 CVT 2 ) — 2 eu> 

Since both matrices are symmetric, the problem is reduced to ensuring that W~ and 
H+ are negative, semi- definite. The eigenvalues of the two matrices are found to be 

Pi, 2 (W - ) = ^ ± \J (C - ) 2 + 16e(/3 2 o»' 2 er 1 2 - 2u>(/3e + 2ao;)r 1 + 2 a; A + e) j 

Pi, 2 (W + ) = ^ ^ ~( + ± ^(C + ) 2 + 16e((J 2 w 2 eT| - 2 uAje + 2 ju> )r 2 - 2 to\ + e)^ , 

where = -2A + 4u)e + Aauri and £+ = 2X + 4cje + 47 ^ 2 . It is evident that negative 
semi-definiteness is ensured if 

/3 2 u; 2 er \ - 2u>(Pe + 2 olw)t\ + 2a>A + e < 0 
6 2 (jj 2 €T$ — 2 u)( 6 e + 2 ju ;)r-2 - 2 u>\ + e < 0 . 

The roots of the two polynomials are 

rf = — ^ + 2 «_ ± 2 yJ k 2 _ + — l/2eo>A^ , 

T± — + 2 *+ ± 2yjK+ + EK+ + l/2£L)\j , 
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where k_ = ua/(3 and k+ -u'y/S. We introduce 

V & = ^6 [£ + 2K “ 2 \/* 2 + e * 

<*> = ^ [e + 2 * + 2 \Z K ' i + £K 

where k — ujafb. Since 

- ■> _l*L i J_ 

“* 6 - 2aw + 4 o> (i) 2 ’ 

for £ < 1 , this ensures > 0 and > 0 . 

Hence, stability is ensured for 

T ~p < ri •< , 

r~ s <t 2 < t + 6 , 

with the solution satisfying 

\j t M % < ~ £ Y1 ■ 0 

3.1.1. Remarks on the Penalty Method for Linear Equations. The results 
stated in Lemma 3.2 may be used to derive the appropriate penalty parameter for 
a large class of linear equations. We consider the general linear ad vect ion- diffusion 
equation , Eq.( 6 ), with the Robin boundary conditions given in Eq.(4)-(5). Solving this 
problem by a penalty method, equivalent to that given by Eq.( 8 ), requires bounds on 
the penalty parameters in order to ensure stability of the scheme. 

In the following, we will give these bounds for reference and will return to the 
numerical validation of these results in Sec. 4. Some of these results may be found in 
[ 3 , 4 , 6 ], but are here given in a more general framework. Note that u) ~ 0(N 2 ), 


- l/2eu\\\\ , 

- l/2eu,\\\\ , 


(i) Hyperbolic Equations (e = 0). 

1. A > 0. Well-posedness is ensured by choosing a > 0 and (3 = 7 = 6 = 
Thus, for this case we will only need bounds on T\. 



r + — 
7 a, 0 — 


00 


The scheme for the hyperbolic case is stable for 


00 > T\ > 


A 

2a )a 
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2 . A < 0 . Well-posedness is ensured by choosing 7 > 0 and a = (3 = S = 0 . 
Thus, for this case we will only need bounds on r 2 . 


r 7,o ~ 


2 o ;7 


<0 = °° 


The scheme for the hyperbolic case is stable for 


OO > T‘2 > 


J*L 

2 u;7 


(ii) Parabolic Equations (A = 0, e > 0). Necessary and sufficient conditions for 
well-posedness may be obtained by choosing the four parameters, a, / 3 , 7 , and 6 y 
properly as stated in Lemma 3.1 [15]. We only state the results for the bounds of T\ y 
since the results for r 2 are equivalent. 

1 . Dirichlet boundary condition, (a > 0, 0 — 0). 




£ 1 
4a (jj 2 


; or,0 


OO 


Stability is ensured for 


OO > T\ > 2 ’ 

4a u) 2 


2. Neumann boundary condition, (a = 0, /? > 0). 


- _ _L + 

T °^ ~ ’ T °^ “ 0U, ’ 


Stability is ensured for 


1 

Tl “ Jkj 

3. Robin boundary condition, (a > 0 , 0 > 0 ). 

= ~g[ £ + 2k - /C 2 + £K] , 

T °'0 = + 2k + 2 >/ k 2 + £*] , 

where k, = cua //?. Stability is ensured for 


T tli ^ ^ T «,/3 • 

(iii) Advection-Diffusion Equations (A / 0, e > 0 ). Again we must ensure well- 
posedness by proper choice of the four parameters, as given in Lemma 3.1. We only 
state the results for the bounds of T\ as the results for r 2 are equivalent. 
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1. Dirichlet boundary condition (a > 0 , 0 = 0 ). 


- - ]*ii J_l_ + 

r “’° - 2a u + 4a a; 2 ’ Tq '° “ 


Stability is ensured for 


| A | 1 el 

°0 >n > + - — 5 

2a a; 4a 

2. Neumann boundary condition (a = 0, /? > 0). 


_ _ 1 _ [2\X& + _ 1 /2jAj^ 

T °-^ _ /3a> V 0 2 ’ T ° ,/3 “ V 0 2 

Stability is ensured for 




~ Tl " (3u 



3. Robin boundary conditions, a > 0, /3 > 0. Results are given in Lemma 3.2. 

3 . 2 . Numerical Tests. As we aim at solving the full non-linear Burgers equa- 
tion, and not the linearized, constant coefficient version, we need to validate the results 
obtained from the linear analysis. We have solved Burgers equation using the scheme 
given by Eq.( 8 ) and employing a standard Legendre collocation method [13, 16]. 

Burgers equation, Eq.(3), has a rightward traveling wave solution (see e.g. [1]) of 
the form 


( 10 ) 


U(Xj t) = — atanh 



+ c 


,x G [- 00 , 00 ] > 0 , 


where the free- stream values 


lim U(x, t) = b — 00 > lim U (x, t) — &oo , 

x — ►— OO r— »-oo 


are associated with the wave-speed, c, and the constant, a > 0 , as 


c = 


^—00 4" ^00 

2 



Following the results in Lemma 3.1 (condition (iv): a = A, /3 = 1, 7 = 0, 6 — 1), we 
expect the non-linear problem to be well-posed for boundary conditions of the type 




dU{-l,t) 

dx 


= 9\{t) 


du(i,t) 

dx 


9i{t) 


where \ > 0 is the value around which we have linearized. In the present study, we 
have used the free-stream value at the inflow, i.e. A = 6 _oo- 
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Since we know an exact solution, the boundary conditions may be given exactly 
at all times using Eq.(10). As initial condition we use 


0) — — a tanh + c . 

The solution is time-stepped using a classical 4th-order Runge-Kutta method, where 
the boundary conditions are imposed at the intermediate time-levels. 

Using the values of the penalty parameters given in Lemma 3.2 results in a stable 
scheme. However, the CFL-number, relating the maximum allowable time step to the 
spatial resolution as 


A, / CFL 

Atniax — i rr i A _1 , A _2 > 

\U As • + £ Ax - 

I I nun 1 min 

will have to be very small in order to ensure stability. Here \U\ signifies the maximum 
absolute value of U. Thus, with the theoretical value of the penalty parameter, the 
proposed method compares unfavorably with the traditional method, due to severe 
time step restrictions. Fortunately, the limits of the penalty parameters, in between 
which asymptotic stability is ensured, are obtained as a result of a conservative energy 
estimate and hence are not very accurate. 

We have used the values of penalty parameter (see Lemma 3.2) as; 



These values are found to lead to a stable scheme, provided eN 2 > 1, In Eq.(3), e 
plays the role of an inverse Reynolds number. The constraint, eN 2 >> 1, simply states 
that increasing the Reynolds number requires increased spatial resolution, which is 
a natural restriction. For advection dominated problems, stability is obtained by 
increasing the penalty parameters towards the values stated in Lemma 3.2. 

With these values of the penalty parameters, we have been able to perform the 
simulations with a CFL number of 4, which is equivalent to what is usually allowed 
when using the traditional method. Thus, by fine-tuning the penalty parameters we 
were able to avoid any effect of the penalty method on the CFL-condition. The follow- 
ing section contains a study of the effect of the penalty method on the C FL - condition 
and guidelines for fine-tuning the penalty parameter for practical applications. 

In Fig. 1 we show the temporal evolution of the traveling wave solution when 
using the proposed scheme as given by Eq.(8). The simulation is done with N = 64 
and e = 0.1. We observe no spurious reflections from the open boundary and the kink 
is seen to travel undisturbed out of the domain. Table 1 shows the error at T — 1.00, 
where the kink has propagated half way through the boundary. It is evident that 
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Fig. 1 . Traveling wave solution of Burgers equation. 

Table 1 

Error in the spectral simulation of Burgers equation using the penalty method. The maximum 
error (Loo) occurs at the boundary. 


N 

L 2 

L oo 

16 

1.07E-02 

3.26E-02 

32 

7.64E-05 

3.50E-04 

64 

3.36E-09 

2.21E-08 

128 

1.56E-11 

7.62E-11 


the proposed scheme maintains the spectral accuracy. The time-step is so small that 
time-stepping errors may be neglected. 

4. CFL-Restrictions for the Penalty Method. As discussed briefly in the 
previous Section, choosing too large a penalty parameter results in severe CFL- 
restrictions. For this reason, it is vital to understand how the penalty method alters the 
eigenvalue spectrum of the operators and consequently changes the CFL- restrict ion. 

In the present section we will study these effects for the linear advection and 
diffusion operators for Legendre collocation methods. For completeness, we will also 
give the results for Chebyshev collocation methods, which are widely used when solving 
non-linear problems. The analysis will consider both 3rd- and 4th-order Runge Kutta 
methods, which are often employed when addressing problems of the type considered 
here. At the end of the section we will compare the results from our linear analysis 
with simulations of the non-linear Burgers equation. 

Consider now the semi-discrete Unear, constant coefficient problem 
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( 11 ) 


(q)t = C N q Xfc 6 n , t > 0 

q = o Xfc e n , t = o , 

= o x t er , t > o 

where q = (q(x 0 ), . - • , q(x/v)) T , k E [0, Cn is the discrete approximation of 

the operator for the interior and Bn determines the appropriate discrete boundary 
conditions. We assume that the semi-discrete approximation is a consistent approxi- 
mation of the continuous problem. A time- differencing scheme, where the boundary 
conditions are enforced exactly at the boundary points, may then be expressed as 

qii+i _ Cn) q n 

B N q n+1 = 0 . 

Here q n signifies the solution vector at time-step n. Thus, for strong stability we must 
require 

\K N (At, C N )\ < 1 • 

However, employing the penalty method changes the time-stepping scheme as 

q n+1 =K N {At,jC N -TB N )q n 

and strong stability is ensured if 

| A”/v(A< , Cn ~ < 1 i 

explaining why the condition depends strongly on the correct choice of the 

penalty-p ar amet er . 

In the following analysis we consider explicit Runge-Kutta time stepping methods, 
which, for time independent operators, may be expressed as 

K^(Ai,CN) = El( At£ Nr , 

t =0 l - 

where p is the order of the scheme. We have for simplicity assumed that the boundary 
conditions are included in the operators. Assuming = SnAjvS^ 1 , where |S;v| and 
1 5^1 are bounded and independent of N , strong stability of the Runge-Kutta schemes 
is obtained if 

m(At, £n)| = Sn j2- ] (AtA N y Sn = < 1 . 

«=0 *• «=0 *• 
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Table 2 


Scaling constants for the advection operator. The proper boundary conditions are of Dirichlet type 


(D)- 


Advection Operator 

Cl 

C c 

A ^ 0 e = 0 

3rd RK 

4th RK 

3rd RK 

4th RK 

D Exact BC 

21 

35 

27 

32 

Penalty BC 

10 

17 

10 

11 


Hence, the problem is reduced to finding the eigenvalue spectrum of the operator £yv 
and choose A t accordingly. 

In the present study we consider the linear advection- diffusion operator; 

r _ x d 9 2 
Ln — A- — b £-^“2 > 
ox ox 1 

with the Robin boundary condition operators 

The boundary conditions for the exact method are enforced through the operator as 
described in [16]. 

In order to compare time-step restrictions as found for the two different ap- 
proaches, we now define the two CFL- like constants, Cl and Cc> as 

At < Cl A Cc 

L ~ XN(N+ l) + eN 2 (N + l) 2 ’ c ~ XN 2 + eN 4 ’ 

where the subscripts refer to Legendre(i) and Chebyshev(C') operators, respectively. 
These constants are determined by solving the eigenvalue problem and calculating the 
maximum At which ensures stability and supplies an upper bound on the time-step. 

Table 2 and 3 shows the calculated values of Cl and Cc for the advection and the 
diffusion operator. The results are the same for the full advection-diffusion operator 
as for the diffusion operator, provided eN 2 > 1, and is therefore omitted. 

It is clear from Table 2 that using the penalty method for enforcing boundary 
conditions on purely advective problems results in a significant reduction of the maxi- 
mum allowable time-step. However, more importantly, Table 3 shows that for problems 
where the diffusion operator dominates the eigenvalue spectrum, the penalty method 
allows for increasing the time-step with as much as 50 %. The effect is most pro- 
nounced when using a 4th-order Runge-Kutta method for time-stepping a Chebyshev 
collocation scheme. 

In order to explain the results in Table 2 and 3, we compare in Fig. 2 the spectrum 
of the Legendre collocation advection (Fig. 2a) and diffusion (Fig. 2b) operators when 
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Table 3 


Scaling constants for the diffusion operator. Results are given for possible combinations of Dirich- 
let (D), Neumann (N) and Robin (R) boundary conditions. 


Diffusion Operator 
A = 0 e > 0 

C 

3rd RK 

L 

4th RK 

C 

3rd RK 

c 

4th RK 

D-D/D-N/D-R 

Exact BC 

99 

109 

53 

58 


Penalty BC 

81 

123 

56 

84 

N-R 

Exact BC 

99 

109 

53 

58 


Penalty BC 

130 

135 

91 

96 

R-R 

Exact BC 

99 

109 

53 

58 


Penalty BC 

130 

141 

93 

97 



• Exaci BC 
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Fig. 2. Eigenvalue spectrum (X — A r + i A t ) for the Legendre advection operator (2a) and the 
Legendre diffusion operator (2b) as obtained by using exact boundary conditions (•) and the penalty 
method (o). 


enforcing Dirichlet boundary conditions through the exact method and the penalty 
method. 

For the advection operator (Fig. 2a) we observe that the effect of the penalty 
method is to introduce an extreme complex conjugate eigenvalue- pair, which domi- 
nates the spectrum and eventually determines the maximum allowable time-step. This 
results in the decreased C F T-number as observed in Table 2. 

The effect on the diffusion operator is more complicated and depends strongly 
on the value of the penalty parameter. As proved by Gottlieb and Lustman [15], 
the diffusion operator with exact Robin boundary conditions has a real, negative and 
distinct eigenvalue spectrum. This property is preserved if a sufficiently large value of 
r is used in the penalty method. However, by decreasing the penalty parameter the 
two most extreme eigenvalues split into two pairs of complex conjugate eigenvalues, 
which move towards the imaginary axis, as r is decreased. In Fig. 2b we show the 
eigenvalue spectrum for the optimal choice of r. The important observation to make 
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is that moduli of these new eigenvalues are smaller than the original extreme negative 
real eigenvalue. Additionally, since the dominating eigenvalue now is complex, it 
clearly becomes advantageous to use the 4th-order Runge-Kutta method due to the 
increased extension of the stability region along ! the imaginary axis as compared to 
The validity of this conclusion is, however, strongly dependent on the proper 
choice of the penalty parameter. The values derived in the previous section do indeed 
ensure asymptotic stability, but with a significant reduction in the maximum allow- 
able CFL-number as a result. Fortunately, as mentioned previously, the limits of the 
penalty parameters are based on a conservative energy estimate and are therefore not 
very accurate. In the following we give the penalty parameters used to obtain the 
results given in Table 2 and 3. These values result in a stable scheme as long as the 
problem is purely advective or eN 2 ;> 1, and allows in most cases for a significant 
increase in the time-step. 


(i) Legendre Collocation Methods 
1. Dirichlet Boundary Conditions. 


JAj 

4 


t = ^N(N + 1) + -^N\N + 1) 2 . 


2. Neumann Boundary Conditions. 

N(N + 1) 


8 


3. Robin Boundary Conditions. 


r — 


Ctfi 


(ii) Chebyshev Collocation Methods 
1. Dirichlet Boundary Conditions. 


r = 1 ^N 2 + — N 4 . 
2 50 


2. Neumann Boundary Conditions. 


N 2 


r = 


3. Robin Boundary Conditions. 


r = with k = 


aN 2 
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Table 4 

Maximum allowable CFL-number obtained from direct numerical simulation of Burgers equation. 



Legendre 

Chebyshev 


3rd RK 

4th RK 

3rd RK 

4th RK 

Exact BC 

3.50 

4.00 

4.25 

4.50 

Penalty BC 

3.00 

4.50 

4.75 

6.75 


Note, that the only difference between the parameter values quoted here and those 
found is Lemma 3.2, is a factor of 1/4 on those terms related to the diffusion operator. 
This reduction is found to lead to optimal time-step restrictions. 

We would like to stress the importance of choosing the appropriate value of the 
penalty parameter. It is our experience, that this is best done by deriving the theoret- 
ical value of this parameter through an analysis similar to that done in Sec. 3.1. This 
leads to a parameter which scales correctly with the resolution and other significant 
parameters. If the time-step restriction is dominated by a viscous time-scale, it is very 
likely that the theoretical estimate leads to severe time-step restrictions. However, 
the theoretical value may often be decreased considerably, and good results may be 
obtained after only a few tests. As we have seen for Burgers equation, decreasing the 
penalty parameter four times leads to acceptable CFL-restrictions. We are not aware 
of any systematic way of determining the optimal factor by which the theoretical value 
should be decreased, but it may usually be determined by trial and error through a 
few tests. 

To conclude our study we have solved Burgers (Eq.(3)) with initial condition 


(12) tf(x,0) = (l-x)(l-x 2 ) , 

and homogeneous Dirichlet boundary conditions. A typical temporal evolution is 
shown in Fig. 3. In Table 4 we show the maximum CFL-number resulting in a stable 
scheme. These results confirm that the results from the linear analysis carries over to 
the non-linear problem. 

5. The Compressible Navier-Stokes Equations. In the present section, we 
obtain energy estimates for the solution to the three-dimensional compressible Navier- 
Stokes equations given in conservation form. Additionally, we derive open boundary 
conditions taking into account the full stress-tensor, and prove well-posedness for the 
continuous problem. The derivations follow the approach introduced in [8, 9]. The 
main difference being that we develop the theory for the conservation form of the 
Navier-Stokes equations and that we include the off-diagonal terms of the stress-tensor 
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Fig. 3. Temporal evolution of Burgers equation with initial conditions given by Eq.(l2). 

in the full derivations. In the second part of this section, we continue by showing how 
to apply the boundary conditions and prove asymptotic stability of the semi-discrete 
scheme. 

Consider now the non-dimensionalized, compressible Navier-Stokes equations given 
in conservation form 

<9q dF dG dU 1 (d dG u dU u \ 

< 13 > a7 + & + ^ + a7 = Ii(^ + -9r + ^rJ ■ 

with x E fi = [— 1, l] 3 . The state vector, q, and the inviscid flux vectors are given as 


p 


pu 


pv 


pw 

pu 


pu 2 + p 


puv 


puw 

pv 

, F = 

puv 

, G = 

pv 2 + p 

, H = 

pvw 

pw 


puw 


pvw 


pw 2 + p 

E 


_ ( E + p)u _ 


. {E + p)v 


. (E + p)w _ 


Here p is the density, u , v, w are the three Cartesian velocity components, E is the 
total energy and p is the pressure. In the remaining part of the paper we will use 
(x,y,z) and (zi,X 2 ,Z 3 ) interchangeably to denote the spatial coordinates. The total 
energy 

£ = P (r+i(u 2 + t , 2 + ™ 2 )) , 

and the pressure is related through the ideal gas law 

P = ( 7 - 1 )pT , 
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where T is the temperature field and 7 = c p /c v is the ratio between the heat capacities 
at constant pressure (c p ) and volume (c^), respectively. 

The viscous flux vectors are given as 


F — 

± V — 


0 




0 

TrX 




Try 

Ty X 


, Gj, 


r yy 

T zx 




r zy 

L r xxU + T yx V + T ZX W + _ 


_ Try^ 4* ^yy^ “1" T zy'W + P7 ( 9 y _ 


r 

0 

- 




r xz 



9 

II 


T yz 





Tzz 




. T XZ U + T yz V + T ZZ W + 



Considering only Newtonian fluids, the stress tensor elements are given as 


TrX 

T yy 

T~zz 


n du 

~ + A 


-- 2fi 


du 


dy 

_ dw f 
2 ^ + A ( 


/ du 


du 


du> A 




/ du 

+ 

du A 

Vd® 

+ 

dy 

+ 

Ih) 

) Try 

— T yx z 

= M 1 

[fry 

dij ’ 

/ du 


du 

+ 

dw\ 



1 

( dw 


duA 

Vd® 

+ 

dy 

dz ) 

) Tyz 

— T zy - 

= H 

<dy 

+ 

d^J ’ 

/ du 


du 


dw\ 




/ dw 


duA 

Vdi 

+ 

dy 

+ 

~dz) 

j Trz 

— Tzx : 

= M 


+ 

dz) 


where p is the dynamic viscosity, A is the bulk viscosity and k is the coefficient of 
thermal conductivity. 

The equations are normalized using the reference values, u re f = uq, p T e f = 
p 0 , Pref = Po'U-oj IVef = Uq/ c v and a reference length X, where (poj^o) is some uni- 
form state, e.g. the ambient free-stream conditions of the flow. This gives a Reynolds 
number as Re = PqUqL/pq and a Prandtl number as Pr = c p po/k 0 . 


5.1. Well-Posedness and Open Boundary Conditions for the Contin- 
uous Problem. Consider the linearized, constant coefficient form of Eq.(13). The 
viscous fluxes are split as 


F„ = F P +F» M + F' M = 


(A + 


L (A + 2/x)u^“ + + f 1W> d^ + Pr J 


dv 


.dw 1 'yk dT 


20 



+ 


0 


0 

*§§ 


\ dw 
* dz 

'‘ft 

+ 

0 

0 


mS 



Auff + M™i . 


G„ = Gp + G x m + G 


M 


M§ 


(A + 2 / 0 g 

L#»«fe+(A + 2#.)i»g + /**S“ + gg 


+ 


0 


0 

Mi 


0 

\ du 

A 3x 

+ 

All 

0 


Mi 

. *»ft + HS . 


. At>§* + /«»£ . 


H, = H p + H i m + H y M = 


M& 

Mi 


(A + 2/x)^ 

[ jitifsf + + (A + 2/i^lf + 



0 


0 


Mi 


0 

+ 

0 

+ 

Mt 


A I 31 

<7X 


Af 3 

ay 




. *»ft + /*»!? . 


Introducing the transformation Jacobians 

dF 


A\ — 


Bn = 


dq 

5Fp 

dqx 


^2 


Bn = 


0G 0H 

^ * M ~~d^ 


dGp 

dq y 


B 33 = 


dH P 

dq z 


_ dG%, 

B,2 ~ 


> ^23 


9G 


M 


, 0HjA o = ( 

dq^ #q v J ’ 13 V 


dl h dll % 

<9q~ dq^ 
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allows for writing Navier-Stokes equations as 


dq 

at 


3 


+ XM» 


dxi 


1 

Re 


3 3 


X] X) 


d 2 q 


It is well known that Navier-Stokes equations, although not of hyperbolic nature, sup- 
port waves very similar to those encountered in the hyperbolic Euler equations. For 
hyperbolic systems, Gottlieb et al. [17] have shown that enforcing the boundary condi- 
tions through the characteristic variables of the system results a stable approximation. 

For Navier-Stokes equations, we linearize around a uniform state, qo, by fixing all 
the matrices. We transform into characteristic variables by diagonalizing A\ through 
a similarity transform A = S~'A\S, where A is the eigenvalue matrix and S and S~ l 
are the matrix of right and left eigenvectors, respectively. These matrices are given in 
the Appendix. Applying this, the symmetrized, linearized set of equations transforms 
into 


(14) 


O 


d 2 K 


2 + ^ A ' dxi Re — ' .. L ,J dx i dxj 

t=l J = l ' J 


where R = S 1 q are the characteristic variables. We have introduced a positive 
definite, symmetrizing diagonal matrix, Q T Q, given as 


O r O = 


1 0 0 0 0 
0 2 0 0 0 
0 0 0 0 
0 0 0 2 0 
0 0 0 0 1 


where cq = y/lPo/po is the uniform state sound speed. Also we define the symmetrized 
matrices 


A* = Q T QS~'A l S , B S %J = Q^QS^BijS . 


The explicit form of the symmetric matrices are given in the Appendix. The charac- 
teristic variables, R = [iZi, R 2 } #3, #4, R$] T > are given as 

r pu - u 0 p + ^ {e + \{ul + vl + wl)p - UQpu - v 0 pv - Wopw'j 1 


R = 


pv - v 0 p 

p - J ^r {e + \p(ul + vl + w'l) - u 0 pu - v 0 pv - w 0 pw^ 


pw - W 0 p 


-(pu - Uop) + {e + \p(ul + Vq + wl) - Uopu - Vo pv - WoP'tvj 
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We are now ready to state the following 


Lemma 5.1. Assume there exists a solution , q ; which is periodic or held at a 
constant value at the y- and z-boundary. If the boundary conditions in the x-direction 
are given such that 


V(t/, z) E Sl v X Sl z 


„ 2 T dK 

J = 1 J 


< 0 


X— — 1 


and the fluid properties are constrained by 


Po > 0 , Aq < 0 , Aq + /io > 0 , 


Pr 


> 0 , 7> 1 , 


then Eq.(14) is a well-posed problem and the solution is bounded as 


1 d_ 

2 dt 




Proof. Construct the energy integral as 


1 d_ 

2 dt 


IIQH.II 2 




dx x Re 


i=i j=i 

3 


/.(-£■ 
f f -\ 

JQ y Jn z 2 

— L f (yy^Ii3 s — ^ dn , 
Rein [hh dx ' ' 3dx >) 


'dx 3 

/ 

1 1 


2 T <9R 
R t a , k ___ Zr t b . 

1 = 1 J 


dydz 


x=-l 


where SI = Sl x X X Sl z . In deriving this expression, we use partial integration 
and assume the solution to be periodic or held at a constant value along the y- and 
z-boundaries, i.e. contributions from these boundaries cancel. This is not a severe 
restriction, as this assumption is valid for a large variety of situations where open 
boundary conditions are applied. 

It is evident that if we can prove 


(15) 



dK r „ s dK\ 

dx t i} dx,) 


dfi > o , 


then well-posedness may be ensured by properly constructing the boundary operator 
at the x-boundary. 
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Since the matrices, Bfj , are all symmetric, Eq.(15) may be rewritten in a block- 
quadratic form as 


where we have introduced 


— !— f R 

2 Re Jq 


T n s R<m > o 







25f, 

B\ 2 

B\3 


'0R 

<m 

0R' 

T 

R = 

.foT’ 


dz . 

, n s = 


2£| 2 

B 23 






^23 

2^3 


We observe that Tt s is a 15 X 15 symmetric matrix, ensuring that the eigenvalue 
spectrum, p(H s ), is real. Hence, if H s is positive semi- definite, Eq.(15) is obeyed. The 
eigenvalue spectrum, p{'H s ) , may be found to be 


p 1 

= Pi = P3 = 0 

P 4 

= 2(^o — Ao) 

PS 

= Pe — 2(Ao + 3/io) 

Pr 

= ps ~ 3/zo — \/ Po + 2(/zo + Ao) 2 

P 9 

= Pio = 3/i 0 + \/^0 + 2 (Mo + Ao) 2 

^>11 

— 7/xo + 4A 0 — \/ Po 4" 4(/io + Ao)(3/zo + 2 Aq) 

P 12 

= 7/io 4- 4Ao \J Po 4(/io 4~ Ao)(3/xo 4- 2Ao) 

P 13 

„ „ ( 2c 2 \2( 7 -1 )*o 

- *>.4 - = 1 (7 _ j )2 + 1 I pi 


Here subscript ’O’ signifies the parameter values in the uniform state around which we 
have linearized. For most real fluids under non-extreme conditions, it is true that 
is positive, Ao is negative and the following relationship is obeyed [18] 


(16) > A 0 + > /zo • 

A simple investigation of the eigenvalues reveals that 7*f 5 is positive semi-definite under 
these conditions. Thus, Eq.(15) is true provided 

IM) > 0 , A 0 < 0 , A o + Mo>0, ^>0, 7 > 1 . 

These conditions are only natural as discussed in [19]. In fact, if they are not obeyed, 
Navier-Stokes equations violates the second law of thermodynamics. 
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We now obtain that well-posedness is ensured under the additional condition 

n l 


V(t/, z) G fit/ X ^ 


and the solution is bounded as 


J = 1 J 


<0 , 


X— — 1 


u 

2 d< 


i^ s -EX(£g^si)" s#i 


where OR = QS 1 q. 


As stated in Lemma 5.1, appropriate boundary conditions at the x-boundary have 
to obey 


~ 2 T dR 

r ’ Ujr -^ ErX ^ 

3 = 1 J 


< 0 . 


We now define 


Q T QG = J2 B h 


3 = 1 


dK 

dxj 


where 


„ _ ^0(7 — 1) d(i Aq + 2//o d (2 ^0 +_/fo / 9R i 2 _ d&A 

Gl “ 2poPr ~dx + 2po dx + po V dy dz ) 

po dRi \o + po d ( 2 
Ct 2 — ^ r 


(17) 


po dx po dy 

kp(j - 1) dCi 

2poCo?x dx 

Po dRt A 0 + po dCi 
G4 — ~ r 


G : 1 — 


G, = 


po dx ' po dz 

^0(7 — l)d£i Aq + 2pp dC 2 Ap + ftp ^ dJ?2 di?4 ^ 


T5 2p 0 Pr dx 2p 0 dx po V dy dz 
where we, for simplicity, have introduced 


Ci — R\ + R$ 


2c 0 


R3 , C2 = R\ - Rs ■ 


7-1 

This allows for rewriting the constraint on the boundary contribution as 

T 1 

O < 0 , 


-io T 
2 ~ 


R t AR- — R r G 
Re 


J-i 


where A is the diagonal eigenvalue matrix obtained from the similarity transform. We 
now reformulate this as 


(18) 


^A- 1 -{eGi) 2 ^ 


<0 , 
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where A; are the wave speeds by which the characteristic variables are advected, as 
given by the diagonal elements of A, and we have introduced e = Re" 1 . This formula- 
tion makes it straightforward to devise inflow-outflow boundary conditions, which are 
maximal dissipative and ensure well-posedness of the complete problem. 

We note in particular that this formulation takes into account the off-diagonal 
terms of the stress tensor, which is neglected in most previous work [8, 9 , 10 ]. These 
terms may be of importance if the artificial boundary is introduced into a strongly 
vortical region of the flow, e.g. a wake flow behind a blunt body. 

Inflow Boundary Conditions. At x = — 1, Eq.( 18 ) becomes 

((|A l | J R 1 - e ^G l ) 2 -(eG,) 2 ) <0 . 

Subsonic Inflow : > 0 , A 2 > 0 , A 3 > 0 , A 4 > 0 , A 5 < 0 

( 19 ) \]R\ — eG\ — 0 

A 2 i? 2 — zG‘2 = 0 

X3R3 — £G 3 = 0 
A 4 fi 4 — £G 4 = 0 

£(j 5 = 0 

Supersonic Inflow : Ai > 0 , A 2 > 0 , A 3 > 0 , A 4 > 0 , A 5 > 0 

(20) \ l R l -eG l =0 

A 2 fi 2 — £G 2 = 0 
A 3 f? 3 — £G 3 = 0 
X4R4 — £G 4 — 0 
A5R5 — eG$ = 0 

Outflow Boundary Conditions. At x = 1, Eq.( 18 ) becomes 

^(|A,|A, - eh^G t ) - (eGi) 2 ^ <0 . 

Subsonic Outflow : Ai > 0 , \ 2 > 0 , A 3 > 0 , A 4 > 0 , A 5 < 0 

(21) eG 2 = 0 

£G 3 = 0 
eG 4 — 0 

I A5 1 7^5 + fGs = 0 
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Supersonic Outflow : Ai > 0 , A 2 > 0 , A 3 > 0 , A 4 > 0 , A 5 > 0 


( 22 ) 


eG\ =0 sG ‘2 — 0 

eG'2 — 0 s G 3 = 0 

° r 

s G% = 0 EG 4 — 0 

eG 4 — 0 zG?) — 0 


We note that for both types of outflow boundary conditions, it is only necessary to 
specify four conditions, since eG 3 = 0 =>• eG\ — — eG$. Due to the special structure 
of G we also observe that adding an extra condition on eG\ does not place extra 
conditions on the solution, since such a condition is redundant. This observation will 
be used later. 

It was shown by Strikwerda [20] that the proper number of boundary conditions 
for an incomplete, parabolic system, like the compressible Navier-Stokes equations, is 
5 in the inflow region and 4 in the outflow region. Our result clearly conforms with 
that. 

We also note that in the limit of infinite Reynolds number, these boundary con- 
ditions converge uniformly toward the well known characteristic boundary conditions 
for the compressible Euler equations [21]. This property is important in order to avoid 
weak boundary layers of the order exp(-xfe) (see [ 8 ]). 

5.2. The Semi-Discrete Scheme. Following the line of thought that led to the 
asymptotically stable scheme for Burgers equation, we propose a Legendre collocation 
scheme for enforcing open boundary conditions to the compressible Navier-Stokes 
equations 


, x dq dF dG dU 1 fdF u , dG u t dU,\ 

(2 % + fa + ^ + “ Re ( dx + dy + dz ) 

-Ti Q-(x)S (n~ (R-S-'giW) - 

-t 2 q+(x)s (n + (r - s-'g 2 (tj) + • 

Here Q~{x) and <3 + (i) are given by Eq.(9) and S is the right eigenvector matrix as 
given in the Appendix. The boundary conditions for the state vector are given through 
the two vectors, gi(t) and g 2 (t), which we for convenience assume to be uniform. The 
four matrices, 7 Z~ , Tl + , Q~ and G + are chosen such as to construct the appropriate 
boundary operator as derived in the previous section. Hence, we have for the inflow 
region 
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where a = 0 for subsonic conditions and a = 1 for supersonic conditions. Likewise we 
define 
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where j3 = 1 for subsonic conditions and f3 = 0 for supersonic outflow conditions. 

We have to choose T\ and 7*2 such that the semi- discrete scheme is asymptotically 
stable. The proper choice is stated in the following Lemma. 


Lemma 5.2. Assume there exists a solution, q, which is periodic or held at a 
constant value at the y- and z-boundary , and that the fluid properties of the uniform 
state, qo, are constrained by 


V- o > 0 , A 0 < 0 , A o + /x o >0, ^ > 0 , 7>1, 

and related as 

~~ > Aq + 2/io A po • 

Pr 

The linearized, constant coefficient version of the scheme given by Eq.(23) is asymp- 
totically stable at the inflow if 


U)K 


^1 + K + VT+^j 


>r x > 



Here 


e k 0 

K ~ — . 

2io Prpo^o 

These results are independent on whether the inflow is subsonic or supersonic. 
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For supersonic outflow 


For subsonic outflow 




i / fi\ i 

- i + v - >r 2 > — . 

U) \ V K / Zu) 


The solution to Eq.(23) is bounded in the form 


Proof Write Eq.(23) is its symmetrized, linearized, constant coefficient version 




dt “ 2 


g2R 

Re ' 3 dxidxj 

1 = 1 J=i J 


-r,g-(i) ^<2 t Q6Tg) 

-r 2 g + (*) (o 7 'Qft+R + ^-Q t Q£+g) , 


where we, without loss of generality, have assumed homogeneous boundary conditions. 
We construct the energy integral, apply the Gauss-Lobatto quadrature rule and partial 
integration to obtain 


(24^|!QR|& =/„ /„ -! 

** 

-r,u J J R T Q T Qn~R- d V dz 

= 1 3 J x=-l 

-Tju [ [ ll 'o 7 Q K 1 R + ^ , 


where we have used the assumption about periodicity or constant value at the y- and 
z-boundary. Additionally, we have introduced e = Re -1 and a>, which is the Legendre 
weight at the endpoints and applied the definition 


QG = 


j-i j 
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Using the Gauss-Lobatto quadrature rule allows for writing 


(25) 



5r7 dn\ 

dxi tJ d Xj j 


dfl = 


/ / (l> fc 

JHy JQ Z y k=\ 


^ 3x, *^x, 


t X—X^ 


+U) 


+ U> 


Y y ^5— /? s ^ 

ax, 13 dx 3 


J x=— 1 


3 3 


vs 

Z.Z. dx . ^ dX} 


dy dz > 0 


J x=l > 


Here x^ signifies the Legendre-Gauss-Lobatto collocation points. The inequality fol- 
lows from the analysis done in the proof for Lemma 5.1, and is ensured provided the 
fluid properties are constrained by 

Mo > 0 > A o <0 , Aq + Mo > 0 , > 0 , 7^1- 

It was shown by Abarbanel and Gottlieb [18] that if a scheme is stable without the 
contributions from the off-diagonal stress-tensor terms, then it will remain so even if 
the these terms are included. This is a consequence of the general relation 


~~ > Ao + 2 mo > Mo j 
Pr 

which roughly gives the relation between the eigenvalues of the normal stress-tensor 
elements and the off-diagonal elements. Thus, it is sufficient to prove stability in the 
absence of the off-diagonal contributions. 

The penalty parameters, T\ and r- 2 , has to be chosen such that the boundary term 
of the energy integral not destroys the stability of the Cauchy- problem. We treat the 
two boundary contributions separately. 


Inflow Condition. The contribution of the boundary term at the inflow (x = — 1) fol- 
lows from combining Eq.(24) and Eq.(25) and neglecting the off-diagonal contributions 
to obtain 

R T (^; -T, u QTm-)R-'R T (l-r,rt-)B t ,g- <0 , 

where 1 is the identity matrix. 

First we note that 

< 9 R t _ OR 


- eu /- 


R s 


dx 2 ^ 22 dx 2 w dx 


dR T n8 OR ^ _ 

£a/ __ ^33^7— < 0 , 


'dx 3 
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since #22 anc ^ ^33 are positive semi-definite with an eigenvalue spectrum given as 


Px(B a 22 ) = = 0 P2 (B‘ 22 ) = P2 (B° 33 ) = f o 

Ps(Bh) = P3{B&) = 2^ p 4 (B* 22 ) = P4 (B° 33 ) = 2^^ . 

Ps(Bh) = (f^ji + l) Ps(B& 3 ) = ((^?F + 0 

Since all matrices are symmetric, the remaining part of the constraint may be expressed 
in block-quadratic form as 

r t h r < 0 , 


where 

r dRl r _ 1 A\ - 2T]U)Q T OR,- -e(l - t x u)B s u 

R_ [ R ’5^J ’ H ~2 — e(l — t\lj)B\ x -2 euB s u 

where we have used Q~ = I. H~ is a 10 X 10 symmetric block-matrix. Similar to 
the approach appHed in Sec. 3.2, we have transformed the problem of stability into 
proving that W", for a suitable value of Ti, is negative semi- definite. The eigenvalue- 
spectrum, p(H~f can be found by doing a LU-decomposition. Since Hr is symmetric, 
the eigenvalues appear as p x (H~) — U l% . 

We will not give the general form of the eigenvalues here, since they are rather 
complicated. However, straightforward but very lengthy algebra shows that all eigen- 
values are negative if T\ is chosen such that 



where 

e k 0 

K — — . 

2 uj Pr/)o^o 

This result is independent of whether the inflow is subsonic or supersonic. 


Outflow Condition. Neglecting the contribution from the off-diagonal terms yields a 
criteria for stability at the outflow (x = 1) 


+t 2 ujQ t Q7Z + ^ R + eR r (1 - t 2 ujQ + ) B s u - 


" dR T 0R 


Similar to the approach followed in the previous part of the proof, we see that the 
contributions from 13 22 an( ^ ^33 are always negative and independently ensure stability. 
We now rewrite the remaining part of the condition at the outflow in block- 


quadratic form; 


R r H+R < 0 , 
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where 


n + 



-A\ - 2 t 2 u>Q t QK + e(l - t 2 u>)I3 s u 
e(l — —2eu)B\ x 


To form we have assumed <? + = I. The- additional boundary condition introduced 
by this replacement is redundant as discussed in Sec. 5.1.2., and, hence, no extra 
restrictions are put on the system by this approach. The eigenvalue spectrum, p(W + ), 
may again be found through a LU-decomposition. We state here only the bounds on 
T ‘2 that ensure negative semi- definiteness of for supersonic outflow 


1 

U) 




For subsonic outflow the bounds become 


1 

CJ 



> T'2 > 


1 

2ct? 


Combining Eq.(24) and Eq.(25), we obtain a bound for the growth of the solution 




dxi 13 dx 3 

t=l J=l J 


dy dz < 0 . □ 


x=x k 


We wish to emphasize that the bounds on Tj and t<i given in Lemma 5.2 remains 
valid in the limit when the Reynolds number approaches infinity. This is easily realized 
by expanding the bounds for e <C 1 to obtain 


OG > T] > 


1 1 

b £ K 

2 u) 8 u 


in the inflow region and 


oo > r 2 > -oo , oo > r 2 > — , 

ZCO 

for supersonic and subsonic outflow, respectively. The linearized, constant coefficient 
version of the Euler equations may be transformed into 5 independent hyperbolic equa- 
tions for which we should expect the bounds on the penalty parameters to be given 
by the results in Sec. 3.1.1. We observe that the bounds given above converge uni- 
formly to the expected values in the limit of vanishing viscosity and, thus, the scheme 
remains stable. The observation that no bounds are necessary on r 2 for supersonic 
outflow simply reflects the fact that no boundary conditions are required for the Euler 
equations at such a boundary. 


32 



5.3. Numerical Tests. The proof given in the previous section is only strictly 
valid for the linearized, constant coefficient version of Navier-Stokes equations. To 
validate the results and show that it carries over to the full non-linear Navier-Stokes 
equations, we have implemented the scheme in an existing spectral code (see [22] for 
details), originally developed for studying two-dimensional compressible flow around 
an infinitely long circular cylinder. 

As spatial approximation scheme was used a standard Fourier- Chebyshev collo- 
cation scheme in polar coordinates, (r, 0), with a 3rd-order Runge-Kutta method for 
time-stepping. 

The new scheme is simple to implement in existing codes, as we only need to apply 
a correction of the flux of the state vector at the boundary. Following the scheme, 
given by Eq.(23), we need to derive the two vectors R and G. The characteristic 
variables are given as 


R\ - ( m T - pu T ) + — , 
co 

r . 2 = p~ L , 

C 0 

R3 = me - pu e , 


R,\ = -(m r - pu r ) + 


co 


where cq is the uniform state sound speed. 
We have for convenience introduced 


u r 


= u 0 ki + v 0 k 2 , uo = u 0 k 2 - v 0 ki , 


which are the radial and azimuthal velocity components, respectively, of the uniform 
state and 


VUr 


= m u k\ + m v ki , m& — Tn u k 2 - m v k\ , 


are the radial and azimuthal components of the momentum of the flow field. Here k = 
(^i> ^ 2 ) signifies an outward pointing normal- vector at the boundary. The linearized 
pressure, p, is given as 


V = (7 - 1) 


E + ~p(uq + Vq) - u 0 m u - Vo m v 


The eigenvalues corresponding to the characteristic functions and determining the 
direction and propagation velocity of the characteristic waves, are 

A! — u r + c 0 , A 2 = A3 = u r , A 4 = u r - c 0 . 
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Following the approach outlined in the previous section, we have likewise derived the 
viscous correction vector, G, at the outer boundary as 


(7 — l)fcp 1 aCi 2 fiodC 2 l/xod-fo 

1 Pr 2po dr 3 po dr 3 po dd ' 

(7 - i)fco 1 dCi 

2 Pr 2c 0 po dr 

„ _ MO j 1 MO ^C'2 

3 po dr 6 po 90 

(7 — l)A:o 1 2 /i 0 

4 ” Pr 2p 0 dr 3 p 0 dr + 3 p 0 90 1 

where again we have defined 

Cl = Jti + R 4 - ~ -R 2 , ( 2 = Rl-R4 ■ 

7-1 

Also, we have J = 1/r, which is the transformation Jacobian from Cartesian to polar 
coordinates. We note that no extra calculation of derivatives is needed in order to 
form the two vectors, since the radial and azimuthal derivatives at the boundary are 
calculated during evaluation of the interior dynamics when employing a global scheme. 
Thus, the only additional requirement is to store values of the derivatives of the state 
vector at the boundary, i.e. the computational requirement for enforcing this new 
method is negligible. 

The boundary conditions are enforced at each intermediate time step of the Runge- 
Kutta method. Simulations were done with a Reynolds number of 100, a Mach number 
of 0.4, the diameter ( D ) of the cylinder being 6.10 cm and the reference temperature 
was 300°A. These parameters ensure that the flow field remains subsonic. The reso- 
lution was 96 Fourier-modes, 72 Chebyshev modes and the radius (i) of the compu- 
tational domain was 20 cylinder diameters. 

As penalty parameters we used 

N 2 2 R x N ' 2 2 

T, = i; i( 1 + *-' /r +^ ) • '> = Ti • 

where N is the number of Chebyshev modes, 2j L is a factor occurring from the radial 
mapping of L into [—1,1] and 

_ eN 2 k 0 
2 Frpouo 

This choice appears naturally from the results stated in Lemma 5.2, and the 
experience gained in Sec. 4, indicating that for dissipative terms we should reduce by 
a factor of 4 in order to obtain the optimal value of r. With this choice of penalty 
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parameters we were able to perform the simulations without any reduction in time- 
step as compared to the exact method of enforcing the boundary conditions. It should 
be mentioned, that in the original code only characteristic boundary conditions for 
the Euler equations were enforced. Comparing with results discussed discussed in Sec. 
4, we observe that for 3rd-order Runge-Kutta we should expect the two methods to 
impose almost equivalent time-step restrictions. This is confirmed by the simulations 
and shows that the results from the simple linear analysis carries over to the full 
non-linear Navier-Stokes equations in this case. 

In Fig. 4 we show contour-plots of the normalized density and the pressure at 
T=143.5, corresponding to approximately 23 shedding cycles. The von Karman vortex 
street is clearly demonstrated, and we observe that the boundary conditions at the 
outflow boundary affect the flow only slightly. The Strouhal number for the shredding 
frequency is found to be St = 0.163, which is in full accordance with experimental 
findings [23] and we observe no spurious frequencies or reflections from the artificial 
boundary back into the flow field (see [22] for a further discussion of this). 


p/po P/Po 



Fig. 4. Contour plots of the normalized density, p/po, cmd the normalized pressure, p/po, at the 
non-dimensional time T=14^.5 for a flow at Re = 100, M = 0.4, D ~ 6.10 cm and 7o = 300°TC. 


6. Concluding Remarks. The purpose of the present paper has been two-fold. 
The first goal has been to develop boundary conditions for wave-dominated problems, 
leading to well-posed total problems. It was argued, that for smooth solutions and the 
kind of operators we have considered here, it is sufficient to consider the problem of 
well-posedness for the linearized, constant coefficient version of the non-linear initial- 
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boundary value problem. Using this allowed for deriving proper boundary conditions to Burgers 
equation and the three-dimensional, compressible Navier-Stokes equations, and these boundary 
conditions were shown to ensure well-posedness of the total problem. It should be stressed that 
the boundary conditions derived for the Navier-Stokes equations takes into account all elements of 
the stress-tensor, and only very light assumptions were made to derive these. Additionally, they 
remain valid even in the limit of vanishing viscosity. 

Having derived appropriate boundary conditions naturally leads to the question of how to 
enforce these in a discrete approximation of the problem. This has been the second, and main, 
contribution of the paper. Recent results [7] on the connection between stability of discrete and 
semi-discrete approximations, suggest that it is sufficient to consider asymptotic stability for the 
semi-discrete approximation. We have only considered Legendre collocation methods here. This 
choice is merely dictated by a wish to obtain analytical results and we have indicated, by numerical 
tests, that all results carry over to Chebyshev collocation operators. The stability proofs for the 
semi-discrete approximations to the linearized, constant coefficient versions of Burgers equation and 
the compressible Navier-Stokes equations are all completed by using the classical energy method. 
We emphasized that the proposed schemes remain stable even in the limit where the problems 
become purely hyperbolic. 

The proposed penalty method changes the eigenvalue spectra of the discrete approximations 
of the operators considerably. In order to understand this, we performed a detailed investigation 
of the effect on the eigenvalue spectra of linear operators. It has been shown that the value of 
the penalty parameter, which is obtained form the theoretical analysis, often implies that the 
maximum allowable time-step compares unfavorable with that allowed through more traditional 
methods. However, we discussed in detail how to remedy this and showed that choosing the 
penalty parameter properly may allow for increasing the maximum time-step with as much as 50%. 
Although we are not aware of a systematic way of determining the optimal value of the penalty 
parameter, we do not see that as any significant disadvantage. Our experience tells that once the 
theoretical values of the penalty parameters are obtained, only a few tests are needed to obtain the 
optimal value. Additionally, this only has to be done once, and since only a few hundred time-steps 
are required to test whether the scheme is stable or not, we consider this an insignificant problem. 

Most of the theoretical results, obtained for linearized, constant coefficient versions of the equa- 
tions, are confirmed by numerical simulations of the full non-linear equations. It is stressed that 
the proposed penalty method is very easy to implement in existing codes, which is an attractive 
feature. 
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Although all results and numerical simulations in this paper are obtained using 

spectral collocation methods, the main conclusions carry over to finite difference/finite 
element methods. The derivation of the proper boundary operators, be that for Burg- 
ers equation or for the compressible, Navier-Stokes equations, is obviously unaffected 
by the choice of spatial approximation method. The proposed penalty method for 
enforcing the boundary conditions may be applied in exactly the same manner as 
discussed here, when using alternative spatial discretization methods. The only dif- 
ference is the value of the penalty parameter, which will depend strongly on the order 
of the method. Thus, applying an other method requires one to derive this penalty 
parameter. This may be done by an approach equivalent to the one utilized here. 
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Appendix: Symmetric Matrices for the Navier-Stokes Equations.. Con- 
sider the linearized, constant coefficient compressible Navier-Stokes equations in con- 
servation form given as 

3 q_. ^ 3 3 

4- ’ 

dt 




dxi Re ~ lJ dxi dx, 

t = l 1 = 1 3=\ J 


The matrix, Ai, diagonalizes under the similarity transform, A = S 1 A\S , where the 
right eigenvector matrix, S, and the left eigenvector matrix, S~ l , are given as 
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Here 

* = f , P = - • 

Introducing this transformation into the Navier-Stokes equations yields 

+ a2R 


dt 1 dx{ Re r—* 13 dx t dx 3 5 

t=i i=i j=i J 
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where R are the characteristic variables and Q T Q is a positive definite, symmetrizing 
diagonal matrix. 

The symmetrized matrices 


At = Q T QS~ l A l S , Btj = Q T QS~'B l3 S 

are given as 


u + c 

0 

0 

0 

0 



V 

c 

0 

0 

0 ' 

0 

2u 

0 

0 

0 



c 

2v 

0 

0 

c 

0 

0 

•7-1 

0 

0 

7 

II 

0 

0 

-fcU 

7-1 

0 

0 

0 

0 

0 

2 u 

0 



0 

0 

0 

2v 

0 

0 

0 

0 

0 

u — c 



_ 0 

c 

0 

0 

V 




w 

0 

0 

c 

0 ' 








0 

2w 

0 

0 

0 







•45 = 

0 

0 ^ArW 

7”1 

0 

0 

> 







c 

0 

0 

2w 

c 








0 

0 

0 

c 

w 






B 


s 

11 — 


2 ~P 


(A + 2 /x) + 9 

0 

2c A 

7-1* 

0 

—(A + 2/z) + 9 

0 

4m 

0 

0 

0 

— Ko 

7-1 

0 

4c 2 Z) 

Tt-Tt 

0 

2c zj 
7-1 * 

0 

0 

0 

4/* 

0 

—(A + 2/x) + 9 

0 

2c >1 
7“1 & 

0 

(A + 2^z) + 9 


ft -\- 9 

0 

2c a 

7-1 9 

0 

—/x + $ 

0 

4(A + 2/x) 

0 

0 

0 

2c A 

n 

4c 2 n 

0 

jc a 


u 


7-1 V 

0 

0 

0 


0 

—fi + 9 

0 

2c 

7-1* 

0 

M + 0 


B 



M + 0 

0 

2c /i 
7-1* 

0 

~M + 0 

0 

4/x 

0 

0 

0 


0 

4c 2 A 

T^TF 0 

0 

1 

^ 1 
I 

^ O 

0 

0 

0 

4(A + 2/x) 

0 

~H + 9 

0 

— ^0 

7-1 

0 

/x + 0 


38 


We have for convenience introduced 
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